Tutorial 2a
Sven Panis
2025-07-06
Tutorial 2a shows various ways to fit Bayesian binary regression models to discrete time-to-event data (section 3). The data are taken from the first experiment of Panis and Schmidt (2016), but using only the no-mask trials (with factor prime = blank, congruent, or incongruent). We performed some prior predictive checks to determine the prior distributions we will employ (section 8). We compare models using WAIC and LOO (section 4), and interpret parameter estimates for one model (section 5). We also plot the logit and cloglog link functions (section 6), and illustrate how various prior distributions for an intercept on the logit and cloglog scales look on the original probability (i.e., hazard) scale (section 7).
1. Load the libraries that we will be using.
pkg <- c("cmdstanr", "standist", "tidyverse", "RColorBrewer", "patchwork",
"brms", "tidybayes", "bayesplot", "future", "parallel", "modelr",
"rstan", "knitr")
lapply(pkg, library, character.only = TRUE)## This is cmdstanr version 0.8.1.9000## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr## - CmdStan path: /Users/spanis/.cmdstan/cmdstan-2.36.0## - CmdStan version: 2.36.0## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: Rcpp
##
## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
##
## Attaching package: 'brms'
##
##
## The following object is masked from 'package:stats':
##
## ar
##
##
##
## Attaching package: 'tidybayes'
##
##
## The following objects are masked from 'package:brms':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t
##
##
## This is bayesplot version 1.11.1
##
## - Online documentation and vignettes at mc-stan.org/bayesplot
##
## - bayesplot theme set to bayesplot::theme_default()
##
## * Does _not_ affect other ggplot2 plots
##
## * See ?bayesplot_theme_set for details on theme setting
##
##
## Attaching package: 'bayesplot'
##
##
## The following object is masked from 'package:brms':
##
## rhat
##
##
## Loading required package: StanHeaders
##
##
## rstan version 2.32.6 (Stan version 2.32.2)
##
##
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
##
##
##
## Attaching package: 'rstan'
##
##
## The following object is masked from 'package:tidyr':
##
## extractSet options.
options(brms.backend = "cmdstanr",
mc.cores = parallel::detectCores(),
future.fork.enable = TRUE,
future.rng.onMisuse = "ignore") ## automatically set in RStudio
supportsMulticore()
detectCores()packageVersion("cmdstanr")
devtools::session_info("rstan")theme settings for ggplot
theme_set(
theme_bw() +
theme(text = element_text(size = 22, face = "bold"),
title = element_text(size = 22, face = "bold"),
legend.position = "bottom")
)
## Set the amount of dodge in figures
pd <- position_dodge(0.7)
pd2 <- position_dodge(1)2. Load and wrangle the person-trial-bin data set that we saved in Tutorial 1a.
ptb_data <- read_csv("Tutorial_1_descriptive_stats/data/inputfile_hazard_modeling.csv")## Rows: 26602 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): condition
## dbl (6): pid, bl, tr, trial, period, event
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(ptb_data)## # A tibble: 6 × 7
## pid bl tr trial condition period event
## <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 1 2 4 4 congruent 1 0
## 2 1 2 4 4 congruent 2 0
## 3 1 2 4 4 congruent 3 0
## 4 1 2 4 4 congruent 4 0
## 5 1 2 4 4 congruent 5 0
## 6 1 2 4 4 congruent 6 0Wrangle the data set.
ptb_data <- ptb_data %>%
# select analysis time range: (200,600] with 10 bins (time bin ranks 6 to 15)
filter(period > 5) %>%
# create categorical predictor for TIME named "timebin" with index coding
mutate(timebin = factor(period,levels=c(6:15)),
# create continuous predictor for TIME named "period_9", centered on bin 9,
period_9 = period - 9,
# create binary variables to indicate each bin
d6 = if_else(period == 6, 1, 0),
d7 = if_else(period == 7, 1, 0),
d8 = if_else(period == 8, 1, 0),
d9 = if_else(period == 9, 1, 0),
d10 = if_else(period == 10, 1, 0),
d11 = if_else(period == 11, 1, 0),
d12 = if_else(period == 12, 1, 0),
d13 = if_else(period == 13, 1, 0),
d14 = if_else(period == 14, 1, 0),
d15 = if_else(period == 15, 1, 0),
# create continuous predictor for trial number named "trial_c", centered on bin 1000 and rescale
trial_c = (trial - 1000)/1000,
# create categorical predictor for trial number named "stage" (early,middle,late) with index coding
stage = ifelse(trial <= 500, 1, ifelse(trial > 1000, 3, 2)),
stage = factor(stage, levels=c(1,2,3)),
# create factor "condition", with "blank" as the reference level
condition = factor(condition, labels = c("blank", "congruent","incongruent")),
# create categorical predictor "prime" with index-coding
prime = ifelse(condition=="blank",1, ifelse(condition=="congruent",2,3)),
prime = factor(prime,levels=c(1,2,3))) %>%
select(pid,event,trial,trial_c,stage,condition,prime,period,period_9,timebin,d6:d15)
head(ptb_data,n=17)## # A tibble: 17 × 20
## pid event trial trial_c stage condition prime period period_9 timebin d6
## <dbl> <dbl> <dbl> <dbl> <fct> <fct> <fct> <dbl> <dbl> <fct> <dbl>
## 1 1 0 4 -0.996 1 congruent 2 6 -3 6 1
## 2 1 0 4 -0.996 1 congruent 2 7 -2 7 0
## 3 1 0 4 -0.996 1 congruent 2 8 -1 8 0
## 4 1 0 4 -0.996 1 congruent 2 9 0 9 0
## 5 1 0 4 -0.996 1 congruent 2 10 1 10 0
## 6 1 0 4 -0.996 1 congruent 2 11 2 11 0
## 7 1 1 4 -0.996 1 congruent 2 12 3 12 0
## 8 1 0 5 -0.995 1 incongru… 3 6 -3 6 1
## 9 1 0 5 -0.995 1 incongru… 3 7 -2 7 0
## 10 1 0 5 -0.995 1 incongru… 3 8 -1 8 0
## 11 1 0 5 -0.995 1 incongru… 3 9 0 9 0
## 12 1 0 5 -0.995 1 incongru… 3 10 1 10 0
## 13 1 0 5 -0.995 1 incongru… 3 11 2 11 0
## 14 1 0 5 -0.995 1 incongru… 3 12 3 12 0
## 15 1 0 5 -0.995 1 incongru… 3 13 4 13 0
## 16 1 0 5 -0.995 1 incongru… 3 14 5 14 0
## 17 1 1 5 -0.995 1 incongru… 3 15 6 15 0
## # ℹ 9 more variables: d7 <dbl>, d8 <dbl>, d9 <dbl>, d10 <dbl>, d11 <dbl>,
## # d12 <dbl>, d13 <dbl>, d14 <dbl>, d15 <dbl>3. Fit hazard models.
Model M0: A null model without experimental predictors
The first model we fit is a “random intercepts” model, where we fit a single grand intercept for each timebin and add random intercepts that vary between participants. This constitutes a general specification of TIME, and can be used if we do not want to make assumptions about how (cloglog-)hazard changes over time (within a trial).
There are two ways to implement such a model. First, we can use the index-coding approach, which provides an intercept for each level of TIME (variable timebin in ptb_data).
Prepare the data file, and specify priors. The skew-normal prior is set for each grand intercept on the cloglog scale, and reflects our prior belief that higher hazard values (above .5) are less likely than lower values (below .5). This prior is plotted in row F of Supplementary Figure 2 (section E in the Supplemental Material). A prior predictive check is performed in section 8.1.
data_M0i <- ptb_data %>% select(pid, event, timebin)
priors_M0i <- c(
set_prior("skew_normal(-1,1,-2)", class = "b"),
set_prior("normal(0, 1)", class = "sd"),
set_prior("lkj(2)", class = "cor")
)Fit model M0i.
plan(multicore)
model_M0i <-
brm(data = data_M0i,
family = bernoulli(link="cloglog"),
formula = event ~ 0 + timebin + (0 + timebin | pid),
prior = priors_M0i,
chains = 4, cores = 4,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
file = "Tutorial_2_Bayesian/models/model_M0i")This took 28 minutes on a MacBook Pro (Sonoma 14.6.1 OS, 18GB Memory, M3 Pro Chip).
model_M0i <- readRDS("Tutorial_2_Bayesian/models/model_M0i.rds")
fixef(model_M0i)## Estimate Est.Error Q2.5 Q97.5
## timebin6 -2.8565057 0.3767531 -3.5094558 -1.9908677
## timebin7 -2.2095218 0.2866088 -2.7448643 -1.5702078
## timebin8 -1.8217410 0.2312341 -2.2779645 -1.3378860
## timebin9 -1.2576980 0.2102598 -1.7099727 -0.8548007
## timebin10 -0.8973223 0.2491912 -1.4494108 -0.4549170
## timebin11 -0.7603867 0.2155859 -1.2527708 -0.3985541
## timebin12 -0.6398048 0.1579879 -0.9997296 -0.3670867
## timebin13 -0.6403525 0.1547766 -0.9892399 -0.3824398
## timebin14 -0.9320306 0.2192195 -1.4376422 -0.5553156
## timebin15 -1.0696255 0.2927871 -1.7113355 -0.5463490#plot(model_M0i)Second, we can fit the same model using a dummy coding approach, as follows.
Prepare the data file, and specify priors.
# select bin indicator variables
data_M0d <- ptb_data %>% select(pid, event, d6:d8, d10:d15)
priors_M0d <- c(
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d6"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d7"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d8"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d10"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d11"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d12"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d13"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d14"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "d15"),
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "Intercept"),
set_prior("normal(0, 1)", class = "sd"),
set_prior("lkj(2)", class = "cor")
)Fit the model with dummy coding (M0d).
plan(multicore)
model_M0d <-
brm(data = data_M0d,
family = bernoulli(link="cloglog"),
formula = event ~ 0 + d6 + d7 + d8 + Intercept + d10 + d11 + d12 + d13 + d14 + d15 +
(d6 + d7 + d8 + 1 + d10 + d11 + d12 + d13 + d14 + d15 | pid),
prior = priors_M0d,
chains = 4, cores = 4,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
file = "Tutorial_2_Bayesian/models/model_M0d")M0d took about 77 minutes.
model_M0d <- readRDS("Tutorial_2_Bayesian/models/model_M0d.rds")
fixef(model_M0d)## Estimate Est.Error Q2.5 Q97.5
## d6 -1.9332829 0.4321629 -2.787111 -1.0573393
## d7 -1.3389718 0.3619834 -2.114674 -0.6658691
## d8 -0.8722236 0.2800040 -1.516462 -0.3925590
## Intercept -1.1067777 0.2591651 -1.664300 -0.6149630
## d10 -0.2136246 0.4512951 -1.255833 0.4375740
## d11 -0.2686736 0.4755249 -1.336390 0.4656902
## d12 -0.1743704 0.4985607 -1.289036 0.5779741
## d13 -0.1899364 0.5138580 -1.299201 0.6065687
## d14 -0.7313939 0.4979868 -1.798154 0.1537556
## d15 -0.6647565 0.4946424 -1.724210 0.1838370Third, we can make assumptions to simplify the model. Instead of using a general specification of TIME (bin rank), we can treat TIME as a continuous variable and make assumptions about how cloglog-hazard changes over time within a trial. For example, we might assume that cloglog-hazard changes in a linear (“period_9”) + quadratic (“period_9_sq”) fashion over time bins within a trial.
Prepare the data file, and specify priors. To find the class b priors, we performed a prior predictive check which is reported in section 8.2.
data_M0c <- ptb_data %>%
select(pid, event, period_9) %>%
mutate(period_9_sq = period_9^2)
priors_M0c <- c(
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "Intercept"),
set_prior("normal(0.3,.25)", class = "b", coef = "period_9"),
set_prior("normal(0,.06)", class= "b", coef="period_9_sq"),
set_prior("normal(0, 1)", class = "sd"),
set_prior("lkj(2)", class = "cor")
)Fit model M0c.
plan(multicore)
model_M0c <-
brm(data = data_M0c,
family = bernoulli(link="cloglog"),
formula = event ~ 0 + Intercept + period_9 + period_9_sq +
(0 + Intercept + period_9 + period_9_sq | pid),
prior = priors_M0c,
chains = 4, cores = 4,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
file = "Tutorial_2_Bayesian/models/model_M0c")Model_M0c took about 19 minutes to run.
model_M0c <- readRDS("Tutorial_2_Bayesian/models/model_M0c.rds")
summary(model_M0c)## Warning: There were 1 divergent transitions after warmup. Increasing
## adapt_delta above 0.999 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: bernoulli
## Links: mu = cloglog
## Formula: event ~ 0 + Intercept + period_9 + period_9_sq + (0 + Intercept + period_9 + period_9_sq | pid)
## Data: data_M0c (Number of observations: 12840)
## Draws: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~pid (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.47 0.19 0.23 0.98 1.00 3410
## sd(period_9) 0.23 0.11 0.12 0.52 1.00 3531
## sd(period_9_sq) 0.05 0.03 0.03 0.12 1.00 3215
## cor(Intercept,period_9) 0.05 0.32 -0.57 0.65 1.00 5826
## cor(Intercept,period_9_sq) -0.31 0.31 -0.81 0.37 1.00 5934
## cor(period_9,period_9_sq) -0.49 0.30 -0.91 0.21 1.00 5149
## Tail_ESS
## sd(Intercept) 4518
## sd(period_9) 4075
## sd(period_9_sq) 4115
## cor(Intercept,period_9) 5288
## cor(Intercept,period_9_sq) 5166
## cor(period_9,period_9_sq) 5129
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.22 0.19 -1.64 -0.85 1.00 3200 3855
## period_9 0.44 0.09 0.23 0.61 1.00 4439 4506
## period_9_sq -0.07 0.02 -0.10 -0.02 1.00 3941 3440
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Model M1. Adding our experimental manipulation.
The next models we fit also include our predictor variable prime type. First, we can use index coding for this categorical predictor.
Prepare the data file, and specify priors.
data_M1i <- ptb_data %>% select(pid, event, timebin, prime)
priors_M1i <- c(
set_prior("skew_normal(-1,1,-2)", class = "b"),
set_prior("normal(0, 1)", class = "sd"),
set_prior("lkj(2)", class = "cor")
)Fit model M1i. The interaction between timebin and prime will specify 30 grand intercepts.
plan(multicore)
model_M1i <-
brm(data = data_M1i,
family = bernoulli(link="cloglog"),
formula = event ~ 0 + timebin:prime +
(0 + timebin:prime | pid),
prior = priors_M1i,
chains = 4, cores = 4,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
file = "Tutorial_2_Bayesian/models/model_M1i")Model M1i took about 124 minutes to fit.
model_M1i <- readRDS("Tutorial_2_Bayesian/models/model_M1i.rds")
fixef(model_M1i)## Estimate Est.Error Q2.5 Q97.5
## timebin6:prime1 -4.4496709 0.7322019 -5.5186672 -2.6172545
## timebin7:prime1 -2.9756681 0.4525940 -3.7424620 -1.9015655
## timebin8:prime1 -1.9072521 0.2783905 -2.4772665 -1.3455780
## timebin9:prime1 -1.1509991 0.2466627 -1.6972177 -0.6858400
## timebin10:prime1 -0.6150386 0.2901719 -1.2785840 -0.1378814
## timebin11:prime1 -0.4684899 0.2254269 -1.0113615 -0.1142685
## timebin12:prime1 -0.4547018 0.2177043 -0.9787083 -0.1061844
## timebin13:prime1 -0.6302033 0.3147734 -1.3489175 -0.1414447
## timebin14:prime1 -0.8417804 0.3360386 -1.6159770 -0.2757566
## timebin15:prime1 -1.2681006 0.4575596 -2.3015065 -0.4789412
## timebin6:prime2 -2.2635497 0.3177752 -2.8664268 -1.5812525
## timebin7:prime2 -1.5347939 0.2431984 -2.0381137 -1.0528518
## timebin8:prime2 -1.2777413 0.2531524 -1.8227535 -0.7926707
## timebin9:prime2 -0.8541505 0.3127625 -1.5475030 -0.3109109
## timebin10:prime2 -0.8885676 0.2999444 -1.5431482 -0.3650001
## timebin11:prime2 -0.8094231 0.2945537 -1.4762308 -0.3187765
## timebin12:prime2 -0.5831442 0.2450133 -1.1480717 -0.1736950
## timebin13:prime2 -0.9003810 0.3698885 -1.7097407 -0.2575006
## timebin14:prime2 -1.2750358 0.4208509 -2.1808997 -0.5112653
## timebin15:prime2 -1.3134203 0.4966061 -2.3831527 -0.4397186
## timebin6:prime3 -2.2520279 0.4140011 -3.0427777 -1.3898838
## timebin7:prime3 -2.0632488 0.4948714 -3.0413115 -1.0857520
## timebin8:prime3 -2.6500772 0.5383075 -3.6544872 -1.5199973
## timebin9:prime3 -2.4139511 0.3149837 -3.0239505 -1.7384690
## timebin10:prime3 -1.8710860 0.4376559 -2.8024945 -1.0426678
## timebin11:prime3 -1.3303723 0.3261964 -2.0448833 -0.7438644
## timebin12:prime3 -1.0256625 0.2455494 -1.5681517 -0.5927057
## timebin13:prime3 -0.9584634 0.2732235 -1.5733905 -0.4806148
## timebin14:prime3 -1.0662756 0.2444208 -1.6111005 -0.6284875
## timebin15:prime3 -1.1332311 0.3770779 -1.9308520 -0.4511871Second, when the shape of the hazard function is rather smooth, as it is for behavioral RT data, one can fit a more parsimonious model by treating TIME as a continuous variable, and use a polynomial specification of the effect of TIME. Thus, if we want to make assumptions about (1) how hazard changes over TIME in the reference condition (blank prime), and (2) how the effect of congruent and incongruent primes change over TIME (relax the proportionality assumption), then we can switch to a dummy coding approach for prime-target congruency (variable “condition”) and treat TIME as a continuous variable (variable “period_9”).
For example, we may assume that hazard can change in a linear + quadratic fashion over time for a blank prime, and that the effects of congruent and incongruent primes relative to blank change in a linear + quadratic fashion, and fit the model called “M1d”.
Get the data, and set priors.
data_M1d <- ptb_data %>%
select(pid, event, period_9, condition) %>%
mutate(period_9_sq = period_9^2)
# check which priors to set for model_M1d
check_priors_M1d <- get_prior(formula = event ~ 0 + Intercept +
condition*period_9 +
condition*period_9_sq +
(1 + condition*period_9 +
condition*period_9_sq | pid),
data = data_M1d,
familiy = bernoulli(link="cloglog"))
priors_M1d <- c(
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "Intercept"),
set_prior("normal(0.3,.25)", class = "b", coef = "period_9"),
set_prior("normal(0,.06)", class= "b", coef = "period_9_sq"),
# A N(0,.4) prior is set for the effects of congruent and incongruent primes
set_prior("normal(0,.4)", class = "b", coef = "conditioncongruent"),
set_prior("normal(0,.4)", class = "b", coef = "conditionincongruent"),
set_prior("normal(0, 1)", class = "sd"),
set_prior("lkj(2)", class = "cor")
)Fit model M1d.
plan(multicore)
model_M1d <-
brm(data = data_M1d,
family = bernoulli(link="cloglog"),
event ~ 0 + Intercept +
condition*period_9 +
condition*period_9_sq +
(1 + condition*period_9 +
condition*period_9_sq | pid),
prior = priors_M1d,
chains = 4, cores = 4,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
file = "Tutorial_2_Bayesian/models/model_M1d")The specification “0 + Intercept + …” removes the default intercept in brm() and adds an explicit Intercept for which we can set the prior ourselves. The variable period_9_sq is a squared version of period_9. Note that duplicate terms in the model formula (e.g., condition) are ignored. Because TIME is centered on bin 9, the Intercept represents the estimated cloglog-hazard in bin 9 for the blank prime condition. Model_M1d took about 184 minutes to run.
model_M1d <- readRDS("Tutorial_2_Bayesian/models/model_M1d.rds")
fixef(model_M1d)## Estimate Est.Error Q2.5
## Intercept -1.15177685 0.25582493 -1.69434975
## conditioncongruent 0.24140379 0.22314135 -0.21795665
## conditionincongruent -0.66334117 0.23005088 -1.04167875
## period_9 0.72996217 0.09864581 0.48435367
## period_9_sq -0.10565127 0.04449052 -0.17370970
## conditioncongruent:period_9 -0.46338867 0.07489813 -0.61617647
## conditionincongruent:period_9 -0.49517690 0.19834295 -0.88920457
## conditioncongruent:period_9_sq 0.08163666 0.02848244 0.02395763
## conditionincongruent:period_9_sq 0.13344507 0.02036826 0.09243871
## Q97.5
## Intercept -0.677335125
## conditioncongruent 0.673433175
## conditionincongruent -0.141682475
## period_9 0.879804825
## period_9_sq 0.001633072
## conditioncongruent:period_9 -0.319271400
## conditionincongruent:period_9 -0.097404735
## conditioncongruent:period_9_sq 0.138767100
## conditionincongruent:period_9_sq 0.173247400One could also fit a model with the variables prime and period_9. However, to include interactions between an index-coded categorical variable and a continuous variable in brm(), one has to switch to the non-linear syntax, as illustrated in the following model formula.
formula = bf(event ~ 0 + a + b * period_9 + c * period_9_sq,
a ~ 0 + prime + (0 + prime |i| pid),
b ~ 0 + prime + (0 + prime |i| pid),
c ~ 0 + prime + (0 + prime |i| pid),
nl = TRUE)The priors could be set as follows.
# check which priors to set
get_prior(formula=bf(event ~ 0 + a + b * period_9 + c * period_9_sq,
a ~ 0 + prime + (0 + prime |i| pid),
b ~ 0 + prime + (0 + prime |i| pid),
c ~ 0 + prime + (0 + prime |i| pid),
nl = TRUE),
data = ptb_data %>%
select(pid, event, period_9, prime) %>%
mutate(period_9_sq = period_9^2),
family = bernoulli(link="cloglog"))## prior class coef group resp dpar nlpar lb ub source
## lkj(1) cor default
## lkj(1) cor pid (vectorized)
## (flat) b a default
## (flat) b prime1 a (vectorized)
## (flat) b prime2 a (vectorized)
## (flat) b prime3 a (vectorized)
## student_t(3, 0, 2.5) sd a 0 default
## student_t(3, 0, 2.5) sd pid a 0 (vectorized)
## student_t(3, 0, 2.5) sd prime1 pid a 0 (vectorized)
## student_t(3, 0, 2.5) sd prime2 pid a 0 (vectorized)
## student_t(3, 0, 2.5) sd prime3 pid a 0 (vectorized)
## (flat) b b default
## (flat) b prime1 b (vectorized)
## (flat) b prime2 b (vectorized)
## (flat) b prime3 b (vectorized)
## student_t(3, 0, 2.5) sd b 0 default
## student_t(3, 0, 2.5) sd pid b 0 (vectorized)
## student_t(3, 0, 2.5) sd prime1 pid b 0 (vectorized)
## student_t(3, 0, 2.5) sd prime2 pid b 0 (vectorized)
## student_t(3, 0, 2.5) sd prime3 pid b 0 (vectorized)
## (flat) b c default
## (flat) b prime1 c (vectorized)
## (flat) b prime2 c (vectorized)
## (flat) b prime3 c (vectorized)
## student_t(3, 0, 2.5) sd c 0 default
## student_t(3, 0, 2.5) sd pid c 0 (vectorized)
## student_t(3, 0, 2.5) sd prime1 pid c 0 (vectorized)
## student_t(3, 0, 2.5) sd prime2 pid c 0 (vectorized)
## student_t(3, 0, 2.5) sd prime3 pid c 0 (vectorized)# set priors
priors <- c(
prior(skew_normal(-1,1,-2), class = b, nlpar = a),
prior(normal(0.3,.25), class = b, nlpar = b),
prior(normal(0,.06), class = b, nlpar = c),
prior(exponential(1), class = sd, group = pid, nlpar = a),
prior(exponential(1), class = sd, group = pid, nlpar = b),
prior(exponential(1), class = sd, group = pid, nlpar = c),
prior(lkj(2), class = cor, group = pid)
)Model M2. Including effects of trial number.
Up until now, we have been working with one time scale, TIME or bin rank within a trial. While many cognitive processes play out on this short time scale (milliseconds to seconds), some play out on longer time scales (minutes to hours), e.g., learning processes.
When we are interested in studying how the hazard of response occurence in our priming experiment also changes on a longer time scale, we can add trial number into the model formula. Here we simply illustrate some possibilities with index and reference coding.
First, we can categorize the predictor trial number by grouping trials in one of a number of “stages” in the experiment (e.g., stage 1 = trials 1 to 500; stage 2 = trials 501 to 1000; stage 3 = trials 1001 and later) and use index coding (variable “stage” in ptb_data).
event ~ 0 + timebin:prime:stage + (0 + timebin:prime:stage | pid)Second, we can treat trial number as a continuous variable, and make assumptions about the way hazard changes with trial number (e.g., linear, quadratic, etc.) for each combination of timebin and prime.
bf(event ~ 0 + a + b * trial_c + c * trial_c_sq,
a ~ 0 + timebin:prime + (0 + timebin:prime |i| pid),
b ~ 0 + timebin:prime + (0 + timebin:prime |i| pid),
c ~ 0 + timebin:prime + (0 + timebin:prime |i| pid),
nl = TRUE)Third, one can use dummy coding for prime type (variable “condition” in ptb_data), and treat TIME (“period_9”) and trial number (“trial_c”) as continous variables.
bf(event ~ 0 + Intercept +
condition*period_9*trial_c +
condition*period_9*trial_c_sq +
condition*period_9_sq*trial_c +
condition*period_9_sq*trial_c_sq +
(1 + condition*period_9*trial_c +
condition*period_9*trial_c^2) +
condition*period_9_sq*trial_c +
condition*period_9_sq*trial_c_sq|pid)4. Compare models using loo and waic.
The predictive accuracy of a set of models can be compared using WAIC and LOO.
model_M0i <- readRDS("Tutorial_2_Bayesian/models/model_M0i.rds")
model_M0d <- readRDS("Tutorial_2_Bayesian/models/model_M0d.rds")
model_M0c <- readRDS("Tutorial_2_Bayesian/models/model_M0c.rds")
model_M1i <- readRDS("Tutorial_2_Bayesian/models/model_M1i.rds")
model_M1d <- readRDS("Tutorial_2_Bayesian/models/model_M1d.rds")Using WAIC and LOO for comparing nonnested models.
model_M0i <- add_criterion(model_M0i, c("loo", "waic"))
model_M0d <- add_criterion(model_M0d, c("loo", "waic"))
model_M0c <- add_criterion(model_M0c, c("loo", "waic"))
model_M1i <- add_criterion(model_M1i, c("loo", "waic"))
model_M1d <- add_criterion(model_M1d, c("loo", "waic"))Compare all five models.
loo_compare(model_M0i, model_M0d, model_M0c, model_M1i, model_M1d, criterion = "loo") %>% print(simplify = F)## elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic
## model_M1i 0.0 0.0 -5111.6 61.7 128.1 3.0 10223.3
## model_M1d -37.4 14.9 -5149.0 62.4 49.5 1.5 10298.0
## model_M0i -420.5 27.3 -5532.1 63.7 49.5 1.0 11064.2
## model_M0d -423.8 27.5 -5535.4 63.8 54.4 1.2 11070.8
## model_M0c -430.7 29.0 -5542.3 63.6 18.3 0.5 11084.7
## se_looic
## model_M1i 123.3
## model_M1d 124.8
## model_M0i 127.5
## model_M0d 127.6
## model_M0c 127.2loo_compare(model_M0i, model_M0d, model_M0c, model_M1i, model_M1d, criterion = "waic") %>% print(simplify = F)## elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic
## model_M1i 0.0 0.0 -5110.5 61.6 126.9 2.9 10221.0
## model_M1d -38.4 14.9 -5148.9 62.4 49.4 1.5 10297.8
## model_M0i -421.5 27.3 -5532.0 63.7 49.4 1.0 11064.0
## model_M0d -424.8 27.4 -5535.3 63.8 54.2 1.2 11070.6
## model_M0c -431.8 28.9 -5542.3 63.6 18.2 0.5 11084.6
## se_waic
## model_M1i 123.3
## model_M1d 124.8
## model_M0i 127.5
## model_M0d 127.6
## model_M0c 127.2model_weights(model_M0i, model_M0d, model_M0c, model_M1i, model_M1d, weights = "loo") %>% round(digits = 3)## model_M0i model_M0d model_M0c model_M1i model_M1d
## 0 0 0 1 0model_weights(model_M0i, model_M0d, model_M0c, model_M1i, model_M1d, weights = "waic") %>% round(digits = 3)## model_M0i model_M0d model_M0c model_M1i model_M1d
## 0 0 0 1 05. Model inference for model M1i.
Diagnostics
The Pareto k estimates can be displayed as follows.
loo(model_M1i)$diagnostics %>%
data.frame() %>%
# attach the `id` values
bind_cols(data_M0i) %>%
mutate(id = 1:n()) %>%
ggplot(aes(x = id, y = pareto_k)) +
geom_point(alpha = 3/4) +
geom_text(data = . %>% filter(pareto_k > .2),
aes(x = id + 2, label = id),
size = 3, hjust = 0) +
theme(panel.grid = element_blank())
Observations with Pareto k values above .7 are problematic.
Posterior distributions of the population-level parameters of M1i
Wrangle the posterior draws.
post <- as_draws_df(model_M1i) %>%
select(-lp__) %>%
as_tibble()
post_summary <- posterior_summary(model_M1i, robust = TRUE)
post_summary[1:30,]## Estimate Est.Error Q2.5 Q97.5
## b_timebin6:prime1 -4.6142800 0.5956642 -5.5186672 -2.6172545
## b_timebin7:prime1 -3.0370900 0.3949646 -3.7424620 -1.9015655
## b_timebin8:prime1 -1.9069150 0.2510412 -2.4772665 -1.3455780
## b_timebin9:prime1 -1.1416050 0.2271491 -1.6972177 -0.6858400
## b_timebin10:prime1 -0.5818410 0.2643394 -1.2785840 -0.1378814
## b_timebin11:prime1 -0.4360185 0.1837638 -1.0113615 -0.1142685
## b_timebin12:prime1 -0.4257980 0.1702447 -0.9787083 -0.1061844
## b_timebin13:prime1 -0.5860980 0.2989470 -1.3489175 -0.1414447
## b_timebin14:prime1 -0.8060520 0.3050790 -1.6159770 -0.2757566
## b_timebin15:prime1 -1.2259300 0.4402358 -2.3015065 -0.4789412
## b_timebin6:prime2 -2.2776150 0.2899002 -2.8664268 -1.5812525
## b_timebin7:prime2 -1.5313350 0.2154885 -2.0381137 -1.0528518
## b_timebin8:prime2 -1.2685700 0.2303812 -1.8227535 -0.7926707
## b_timebin9:prime2 -0.8250290 0.2849350 -1.5475030 -0.3109109
## b_timebin10:prime2 -0.8661925 0.2812425 -1.5431482 -0.3650001
## b_timebin11:prime2 -0.7774375 0.2725523 -1.4762308 -0.3187765
## b_timebin12:prime2 -0.5605580 0.2181727 -1.1480717 -0.1736950
## b_timebin13:prime2 -0.8646980 0.3544140 -1.7097407 -0.2575006
## b_timebin14:prime2 -1.2479250 0.4003391 -2.1808997 -0.5112653
## b_timebin15:prime2 -1.2724250 0.4825048 -2.3831527 -0.4397186
## b_timebin6:prime3 -2.2622500 0.3914657 -3.0427777 -1.3898838
## b_timebin7:prime3 -2.0636400 0.4873454 -3.0413115 -1.0857520
## b_timebin8:prime3 -2.6676000 0.5286359 -3.6544872 -1.5199973
## b_timebin9:prime3 -2.4215700 0.2800557 -3.0239505 -1.7384690
## b_timebin10:prime3 -1.8555650 0.4168923 -2.8024945 -1.0426678
## b_timebin11:prime3 -1.3114950 0.2997150 -2.0448833 -0.7438644
## b_timebin12:prime3 -1.0061350 0.2263708 -1.5681517 -0.5927057
## b_timebin13:prime3 -0.9334835 0.2431501 -1.5733905 -0.4806148
## b_timebin14:prime3 -1.0480200 0.2188021 -1.6111005 -0.6284875
## b_timebin15:prime3 -1.1136550 0.3581458 -1.9308520 -0.4511871post_qi_b <- post %>%
select(starts_with("b_")) %>%
pivot_longer(everything()) %>%
group_by(name) %>%
median_qi(value) %>%
arrange(name)
head(post_qi_b) # 30 "fixed" effects ## # A tibble: 6 × 7
## name value .lower .upper .width .point .interval
## <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 b_timebin10:prime1 -0.582 -1.28 -0.138 0.95 median qi
## 2 b_timebin10:prime2 -0.866 -1.54 -0.365 0.95 median qi
## 3 b_timebin10:prime3 -1.86 -2.80 -1.04 0.95 median qi
## 4 b_timebin11:prime1 -0.436 -1.01 -0.114 0.95 median qi
## 5 b_timebin11:prime2 -0.777 -1.48 -0.319 0.95 median qi
## 6 b_timebin11:prime3 -1.31 -2.04 -0.744 0.95 median qiVisualise fixed effects.
tidy_fixed <- post %>%
select(starts_with("b_"), .chain, .iteration, .draw) %>%
rename(chain=.chain, iter=.iteration, draw=.draw) %>%
pivot_longer(-c(chain, draw, iter)) %>%
mutate(timebin = str_sub(name,10,11),
timebin = factor(str_remove(timebin,":"),levels=c(6:15)),
condition = str_sub(name,17,18),
condition = factor(str_remove(condition,"e"),
levels=c(1,2,3),
labels=c("blank","congruent","incongruent")))
head(tidy_fixed)## # A tibble: 6 × 7
## chain iter draw name value timebin condition
## <int> <int> <int> <chr> <dbl> <fct> <fct>
## 1 1 1 1 b_timebin6:prime1 -4.63 6 blank
## 2 1 1 1 b_timebin7:prime1 -2.75 7 blank
## 3 1 1 1 b_timebin8:prime1 -1.62 8 blank
## 4 1 1 1 b_timebin9:prime1 -0.900 9 blank
## 5 1 1 1 b_timebin10:prime1 -0.151 10 blank
## 6 1 1 1 b_timebin11:prime1 -0.449 11 blanktail(tidy_fixed)## # A tibble: 6 × 7
## chain iter draw name value timebin condition
## <int> <int> <int> <chr> <dbl> <fct> <fct>
## 1 4 2000 8000 b_timebin10:prime3 -2.22 10 incongruent
## 2 4 2000 8000 b_timebin11:prime3 -0.981 11 incongruent
## 3 4 2000 8000 b_timebin12:prime3 -1.10 12 incongruent
## 4 4 2000 8000 b_timebin13:prime3 -1.05 13 incongruent
## 5 4 2000 8000 b_timebin14:prime3 -1.07 14 incongruent
## 6 4 2000 8000 b_timebin15:prime3 -1.19 15 incongruent# plot
p_tidy_fixed <- ggplot(tidy_fixed, aes(x = timebin, y = value, fill=condition)) +
stat_halfeye(point_interval="median_qi",
.width = c(0.80,0.95),
alpha=0.7) +
labs(title = 'Posterior distributions for population-level\neffects in Model M1i',
x = "time bin", y = "cloglog-hazard") +
scale_fill_brewer(palette = "Dark2") +
scale_x_discrete(labels = str_c("(",(c(6:15)-1)*40,",",c(6:15)*40,"]"), breaks = 6:15) +
theme(axis.text.x = element_text(angle=90)) +
facet_wrap(~condition)
p_tidy_fixed
Save the posterior distribution plot.
ggsave("Tutorial_2_Bayesian/figures/M1i_postdistr.png", width = 10, height = 8, dpi = 800)Expected value of the posterior predictive distribution
Following Heiss (2021), we can plot the expected value of the posterior predictive distribution – the predicted hazard values – for the “average” participant at the population level, and for each participant in the data set.
First, at the population level, using add_epred_draws().
dat_M1i <- as_tibble(model_M1i$data)
epreds_grand <- dat_M1i %>%
data_grid(timebin, prime) %>%
add_epred_draws(model_M1i,
re_formula = NA) %>% # ignore random effects
mutate(prime = factor(prime,
levels=c(1,2,3),
labels=c("blank","congruent","incongruent"))) %>%
ungroup()Summarize and plot predicted hazard values.
epreds_grand %>%
group_by(timebin,prime) %>%
median_qi(.width=0.95)## # A tibble: 30 × 8
## timebin prime .epred .lower .upper .width .point .interval
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 blank 0.00986 0.00400 0.0704 0.95 median qi
## 2 6 congruent 0.0974 0.0553 0.186 0.95 median qi
## 3 6 incongruent 0.0989 0.0466 0.221 0.95 median qi
## 4 7 blank 0.0468 0.0234 0.139 0.95 median qi
## 5 7 congruent 0.194 0.122 0.295 0.95 median qi
## 6 7 incongruent 0.119 0.0466 0.287 0.95 median qi
## 7 8 blank 0.138 0.0805 0.229 0.95 median qi
## 8 8 congruent 0.245 0.149 0.364 0.95 median qi
## 9 8 incongruent 0.0671 0.0255 0.196 0.95 median qi
## 10 9 blank 0.273 0.167 0.396 0.95 median qi
## # ℹ 20 more rows p1 <- ggplot(epreds_grand, aes(x=timebin, y=.epred,
fill=prime, color=prime)) +
stat_lineribbon(point_interval="median_qi",.width=c(.8,.95),alpha = 0.5) +
scale_fill_brewer(palette = "Dark2") +
scale_color_brewer(palette = "Dark2") +
scale_x_discrete(labels = str_c("(",(c(6:15)-1)*40,",",c(6:15)*40,"]"), breaks = 6:15) +
theme(axis.text.x = element_text(angle=90)) +
scale_y_continuous(limits=c(0,.7)) +
labs(x = "time bin", y = "predicted hazard")
p1## Warning: Removed 19 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).
Again, for the “average” participant at the population level, but now using posterior draws.
tidy_fixed %>%
mutate(haz = 1-exp((-1)*exp(value))) %>% # inverse cloglog
ggplot(aes(x=timebin, y=haz,
fill=condition, color=condition)) +
stat_lineribbon(point_interval="median_qi",
.width=.95,
alpha = 0.5) +
scale_fill_brewer(palette = "Dark2") +
scale_color_brewer(palette = "Dark2") +
scale_x_discrete(labels = str_c("(",(c(6:15)-1)*40,",",c(6:15)*40,"]"), breaks = 6:15) +
theme(axis.text.x = element_text(angle=90)) +
scale_y_continuous(limits=c(0,.6)) +
labs(x = "time bin", y = "predicted hazard") ## Warning: Removed 743 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).
Second, for each participant in the data set.
epreds_pid <- dat_M1i %>%
data_grid(pid,timebin, prime) %>%
add_epred_draws(model_M1i,
re_formula = NULL) %>% # include random effects
mutate(prime = factor(prime,
levels=c(1,2,3),
labels=c("blank","congruent","incongruent"))) %>%
ungroup()Summarize and plot predicted hazard values.
epreds_pid %>%
group_by(pid,timebin,prime) %>%
median_qi(.width=0.95)## # A tibble: 180 × 9
## pid timebin prime .epred .lower .upper .width .point .interval
## <dbl> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 1 6 blank 0.00948 0.00354 0.0271 0.95 median qi
## 2 1 6 congruent 0.0744 0.0400 0.123 0.95 median qi
## 3 1 6 incongruent 0.0742 0.0393 0.125 0.95 median qi
## 4 1 7 blank 0.0665 0.0390 0.107 0.95 median qi
## 5 1 7 congruent 0.176 0.117 0.246 0.95 median qi
## 6 1 7 incongruent 0.101 0.0557 0.166 0.95 median qi
## 7 1 8 blank 0.222 0.170 0.284 0.95 median qi
## 8 1 8 congruent 0.272 0.194 0.363 0.95 median qi
## 9 1 8 incongruent 0.0375 0.0137 0.0812 0.95 median qi
## 10 1 9 blank 0.404 0.331 0.481 0.95 median qi
## # ℹ 170 more rowsp2 <- ggplot(epreds_pid, aes(x=timebin, y=.epred,
fill=prime, color=prime)) +
stat_lineribbon(point_interval="median_qi",
.width=c(.8,.95),
alpha = 0.5) +
scale_fill_brewer(palette = "Dark2") +
scale_color_brewer(palette = "Dark2") +
scale_x_discrete(labels = str_c("(",(c(6:15)-1)*40,",",c(6:15)*40,"]"), breaks = 6:15) +
theme(axis.text.x = element_text(angle=90)) +
scale_y_continuous(limits=c(0,1)) +
labs(x = "time bin", y = "predicted hazard") +
# ggtitle("Each participant (N = 6)") +
facet_wrap(~pid, ncol=3)
p2
Combine the plots.
# combine 2 plots
p1_theme <- p1 + theme(legend.position = "none")
(p1_theme / p2) +
plot_annotation(tag_levels = 'A') +
plot_layout(guides = "collect",
axes = "collect_x") &
theme(legend.position = "bottom")Save the combined plot.
# save combined plot
ggsave("Tutorial_2_Bayesian/figures/M1i_pred_combined.png", width = 10, height = 12, dpi = 800)Average marginal effects (AMEs)
Grand AMEs
Following Heiss (2021), we are actually interested in the difference in predicted hazard between congruent and blank primes on the one hand, and between incongruent and blank primes on the other hand, for each time bin. When based on the grand mean, this is known as the grand average marginal effect (AME).
epreds_grand_diffs <- epreds_grand %>%
pivot_wider(id_cols = c(timebin, .draw),
names_from = "prime",
values_from = ".epred") %>% #80000 x 5
mutate(`congruent minus blank` = congruent - blank,
`incongruent minus blank` = incongruent - blank) %>%
select(-c(blank,congruent,incongruent)) %>%
pivot_longer(cols = c(`congruent minus blank`,`incongruent minus blank`),
names_to = "contrast",
values_to = "diff")
epreds_grand_diffs # 160 000 rows ## # A tibble: 160,000 × 4
## timebin .draw contrast diff
## <fct> <int> <chr> <dbl>
## 1 6 1 congruent minus blank 0.0834
## 2 6 1 incongruent minus blank 0.0805
## 3 6 2 congruent minus blank 0.149
## 4 6 2 incongruent minus blank 0.0342
## 5 6 3 congruent minus blank 0.0982
## 6 6 3 incongruent minus blank 0.0861
## 7 6 4 congruent minus blank 0.0592
## 8 6 4 incongruent minus blank 0.0476
## 9 6 5 congruent minus blank 0.0403
## 10 6 5 incongruent minus blank 0.0709
## # ℹ 159,990 more rowsSummarize.
table <- epreds_grand_diffs %>%
group_by(contrast,timebin) %>%
mean_qi(.width = c(.95)) %>%
arrange(contrast,timebin,.width) %>%
select(-c(.width,.point,.interval))Save summary as table.
write_csv(table, file="Tutorial_2_Bayesian/tables/grand_AMEs.csv")Show table.
kable(table, caption = "Point (mean) and 95% credible interval summary of estimated differences in hazard, for each time bin and contrast, in the average participant.")## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")## Warning in attr(x, "format"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")| contrast | timebin | diff | .lower | .upper |
|---|---|---|---|---|
| congruent minus blank | 6 | 0.0869377 | 0.0165561 | 0.1703592 |
| congruent minus blank | 7 | 0.1427225 | 0.0378084 | 0.2465272 |
| congruent minus blank | 8 | 0.1057670 | -0.0213992 | 0.2350213 |
| congruent minus blank | 9 | 0.0783076 | -0.1156226 | 0.2712171 |
| congruent minus blank | 10 | -0.0794645 | -0.2996066 | 0.1516519 |
| congruent minus blank | 11 | -0.1026104 | -0.3077325 | 0.1108873 |
| congruent minus blank | 12 | -0.0409448 | -0.2413813 | 0.1651772 |
| congruent minus blank | 13 | -0.0755303 | -0.3256924 | 0.1959752 |
| congruent minus blank | 14 | -0.1021703 | -0.3335525 | 0.1549900 |
| congruent minus blank | 15 | -0.0070710 | -0.2658903 | 0.2655349 |
| incongruent minus blank | 6 | 0.0911091 | 0.0146335 | 0.2056552 |
| incongruent minus blank | 7 | 0.0761513 | -0.0349575 | 0.2333228 |
| incongruent minus blank | 8 | -0.0645286 | -0.1692114 | 0.0708865 |
| incongruent minus blank | 9 | -0.1864419 | -0.3117364 | -0.0646903 |
| incongruent minus blank | 10 | -0.2704697 | -0.4530995 | -0.0430156 |
| incongruent minus blank | 11 | -0.2285755 | -0.3958541 | -0.0253659 |
| incongruent minus blank | 12 | -0.1668673 | -0.3309868 | 0.0232611 |
| incongruent minus blank | 13 | -0.0956633 | -0.3147249 | 0.1363460 |
| incongruent minus blank | 14 | -0.0628414 | -0.2734297 | 0.1462738 |
| incongruent minus blank | 15 | 0.0257797 | -0.2258818 | 0.2696414 |
Plot.
p11 <- ggplot(epreds_grand_diffs, aes(x=timebin, y=diff,
fill=contrast)) +
stat_lineribbon(alpha = 0.4,
point_interval = "mean_qi",
.width=c(.8,.95)) +
geom_hline(yintercept = 0,
color = "red",
lty = 3) +
scale_y_continuous(limits=c(-0.5,0.3)) +
scale_x_discrete(labels = str_c("(",(c(6:15)-1)*40,",",c(6:15)*40,"]"), breaks = 6:15) +
scale_fill_brewer(palette = "Dark2") +
scale_color_brewer(palette = "Dark2") +
theme(axis.text.x = element_text(angle=90)) +
labs(x = "time bin", y = "difference in predicted hazard") +
#ggtitle("Grand average marginal effect") +
facet_wrap(~contrast)
p11## Warning: Removed 357 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).## Warning: Removed 262 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).
Note that these grand AMEs can also be calculated based on the posterior draws.
post %>% # 8000 x 30 = 240000
mutate_all(function(x){1-exp((-1)*exp(x))}) %>%
mutate(bin6_cb = `b_timebin6:prime2` - `b_timebin6:prime1`,
bin7_cb = `b_timebin7:prime2` - `b_timebin7:prime1`,
bin8_cb = `b_timebin8:prime2` - `b_timebin8:prime1`,
bin9_cb = `b_timebin9:prime2` - `b_timebin9:prime1`,
bin10_cb = `b_timebin10:prime2` - `b_timebin10:prime1`,
bin11_cb = `b_timebin11:prime2` - `b_timebin11:prime1`,
bin12_cb = `b_timebin12:prime2` - `b_timebin12:prime1`,
bin13_cb = `b_timebin13:prime2` - `b_timebin13:prime1`,
bin14_cb = `b_timebin14:prime2` - `b_timebin14:prime1`,
bin15_cb = `b_timebin15:prime2` - `b_timebin15:prime1`,
bin6_ib = `b_timebin6:prime3` - `b_timebin6:prime1`,
bin7_ib = `b_timebin7:prime3` - `b_timebin7:prime1`,
bin8_ib = `b_timebin8:prime3` - `b_timebin8:prime1`,
bin9_ib = `b_timebin9:prime3` - `b_timebin9:prime1`,
bin10_ib = `b_timebin10:prime3` - `b_timebin10:prime1`,
bin11_ib = `b_timebin11:prime3` - `b_timebin11:prime1`,
bin12_ib = `b_timebin12:prime3` - `b_timebin12:prime1`,
bin13_ib = `b_timebin13:prime3` - `b_timebin13:prime1`,
bin14_ib = `b_timebin14:prime3` - `b_timebin14:prime1`,
bin15_ib = `b_timebin15:prime3` - `b_timebin15:prime1`) %>%
select(starts_with("bin")) %>% # 8000 x 20
pivot_longer(cols = bin6_cb:bin15_ib,
names_to = "condition",
values_to = "diff") %>%
mutate(bin = str_sub(condition,4,5),
bin = str_remove(bin,"_"),
bin = factor(bin, levels=c(6:15)),
comp = str_sub(condition,6,8),
comp = str_remove(comp,"_"),
contrast = ifelse(comp == "cb",
"congruent minus blank", "incongruent minus blank")) %>%
ggplot(aes(x=bin, y=diff, fill=contrast)) +
stat_lineribbon(alpha = 0.5,
point_interval = "mean_qi",
.width=c(.8,.95)) +
geom_hline(yintercept = 0,
color = "red",
lty = 3) +
scale_y_continuous(limits=c(-0.5,0.3)) +
scale_x_discrete(labels = str_c("(",(c(6:15)-1)*40,",",c(6:15)*40,"]"), breaks = 6:15) +
scale_fill_brewer(palette = "Dark2") +
scale_color_brewer(palette = "Dark2") +
theme(axis.text.x = element_text(angle=90)) +
labs(y = "difference in predicted hazard", x = "timebin") +
ggtitle("Grand average marginal effect") +
facet_wrap(~contrast)## Warning: Removed 357 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).## Warning: Removed 262 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).
Subject-specific average marginal effects
To calculate the subject-specific AMEs for each time bin, we create contrasts between the conditional means.
epreds_pid_diffs <- epreds_pid %>%
pivot_wider(id_cols = c(pid,timebin, .draw),
names_from = "prime",
values_from = ".epred") %>% #80000 x 5
mutate(`congruent minus blank` = congruent - blank,
`incongruent minus blank` = incongruent - blank) %>%
select(-c(blank,congruent,incongruent)) %>%
pivot_longer(cols = c(`congruent minus blank`,`incongruent minus blank`),
names_to = "contrast",
values_to = "diff")
epreds_pid_diffs # 160 000 rows ## # A tibble: 960,000 × 5
## pid timebin .draw contrast diff
## <dbl> <fct> <int> <chr> <dbl>
## 1 1 6 1 congruent minus blank 0.0567
## 2 1 6 1 incongruent minus blank 0.117
## 3 1 6 2 congruent minus blank 0.0370
## 4 1 6 2 incongruent minus blank 0.0111
## 5 1 6 3 congruent minus blank 0.0603
## 6 1 6 3 incongruent minus blank 0.0643
## 7 1 6 4 congruent minus blank 0.113
## 8 1 6 4 incongruent minus blank 0.0344
## 9 1 6 5 congruent minus blank 0.0286
## 10 1 6 5 incongruent minus blank 0.0748
## # ℹ 959,990 more rowsSummarize the contrasts in predicted hazard values.
table <- epreds_pid_diffs %>%
group_by(pid,timebin,contrast) %>%
mean_qi(.width = .95) %>%
arrange(contrast,timebin,.width) %>%
select(-c(.width,.point,.interval))Save as table.
write_csv(table, file="Tutorial_2_Bayesian/tables/pid_AMEs.csv")Show table.
kable(table, caption = "Point (mean) and 95% credible interval summary of estimated differences in hazard, for each time bin and contrast, in the average participant.")## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")## Warning in attr(x, "format"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")| pid | timebin | contrast | diff | .lower | .upper |
|---|---|---|---|---|---|
| 1 | 6 | congruent minus blank | 0.0651522 | 0.0262784 | 0.1123884 |
| 2 | 6 | congruent minus blank | 0.0525837 | 0.0190572 | 0.0956194 |
| 3 | 6 | congruent minus blank | 0.1198058 | 0.0697290 | 0.1852845 |
| 4 | 6 | congruent minus blank | 0.0560210 | 0.0221626 | 0.1006466 |
| 5 | 6 | congruent minus blank | 0.1661312 | 0.0978048 | 0.2491776 |
| 6 | 6 | congruent minus blank | 0.0660552 | 0.0283218 | 0.1135927 |
| 1 | 7 | congruent minus blank | 0.1097783 | 0.0390644 | 0.1842512 |
| 2 | 7 | congruent minus blank | 0.1234487 | 0.0660270 | 0.1892208 |
| 3 | 7 | congruent minus blank | 0.1544407 | 0.0847781 | 0.2303358 |
| 4 | 7 | congruent minus blank | 0.1069777 | 0.0424673 | 0.1754384 |
| 5 | 7 | congruent minus blank | 0.2919527 | 0.1814899 | 0.4088326 |
| 6 | 7 | congruent minus blank | 0.1444816 | 0.0790017 | 0.2170399 |
| 1 | 8 | congruent minus blank | 0.0509000 | -0.0502398 | 0.1539529 |
| 2 | 8 | congruent minus blank | 0.0327057 | -0.0376418 | 0.1123674 |
| 3 | 8 | congruent minus blank | 0.0849463 | -0.0064113 | 0.1827095 |
| 4 | 8 | congruent minus blank | 0.0111541 | -0.0828538 | 0.1131906 |
| 5 | 8 | congruent minus blank | 0.1743302 | 0.0566536 | 0.3054037 |
| 6 | 8 | congruent minus blank | 0.2908410 | 0.1906566 | 0.3943589 |
| 1 | 9 | congruent minus blank | 0.0174080 | -0.1177691 | 0.1540204 |
| 2 | 9 | congruent minus blank | -0.1678800 | -0.2636435 | -0.0678977 |
| 3 | 9 | congruent minus blank | 0.1520304 | 0.0192695 | 0.2882150 |
| 4 | 9 | congruent minus blank | 0.0425950 | -0.0929289 | 0.1769888 |
| 5 | 9 | congruent minus blank | 0.2233430 | 0.0583618 | 0.3903858 |
| 6 | 9 | congruent minus blank | 0.2850183 | 0.1564384 | 0.4186883 |
| 1 | 10 | congruent minus blank | -0.2995293 | -0.4685392 | -0.1277672 |
| 2 | 10 | congruent minus blank | -0.2459417 | -0.3762226 | -0.1159232 |
| 3 | 10 | congruent minus blank | -0.1980731 | -0.3543228 | -0.0300314 |
| 4 | 10 | congruent minus blank | 0.0417529 | -0.1776587 | 0.2287208 |
| 5 | 10 | congruent minus blank | -0.0916297 | -0.2845071 | 0.1208433 |
| 6 | 10 | congruent minus blank | 0.1167459 | -0.0413109 | 0.2876520 |
| 1 | 11 | congruent minus blank | -0.1387265 | -0.3847258 | 0.1038641 |
| 2 | 11 | congruent minus blank | 0.0197533 | -0.1729061 | 0.2078299 |
| 3 | 11 | congruent minus blank | -0.1437507 | -0.3341952 | 0.0462079 |
| 4 | 11 | congruent minus blank | -0.0902666 | -0.3944417 | 0.2407306 |
| 5 | 11 | congruent minus blank | -0.1483164 | -0.3740691 | 0.1020338 |
| 6 | 11 | congruent minus blank | -0.0772787 | -0.2618849 | 0.1180696 |
| 1 | 12 | congruent minus blank | -0.0543442 | -0.3113425 | 0.2063196 |
| 2 | 12 | congruent minus blank | -0.1016176 | -0.3410021 | 0.1246328 |
| 3 | 12 | congruent minus blank | -0.0326060 | -0.2276414 | 0.1821579 |
| 4 | 12 | congruent minus blank | -0.0347405 | -0.3260749 | 0.2593652 |
| 5 | 12 | congruent minus blank | -0.0709364 | -0.3034481 | 0.1857469 |
| 6 | 12 | congruent minus blank | -0.0211521 | -0.2245901 | 0.1717465 |
| 1 | 13 | congruent minus blank | 0.1341784 | -0.2498373 | 0.5733570 |
| 2 | 13 | congruent minus blank | 0.1062506 | -0.2485932 | 0.4830828 |
| 3 | 13 | congruent minus blank | -0.1348162 | -0.4228327 | 0.1650868 |
| 4 | 13 | congruent minus blank | -0.0013988 | -0.4101252 | 0.4485768 |
| 5 | 13 | congruent minus blank | -0.0830838 | -0.4215040 | 0.3231655 |
| 6 | 13 | congruent minus blank | -0.2705035 | -0.5274898 | 0.0092309 |
| 1 | 14 | congruent minus blank | -0.0331333 | -0.4221919 | 0.5133149 |
| 2 | 14 | congruent minus blank | -0.0242911 | -0.4152219 | 0.4636386 |
| 3 | 14 | congruent minus blank | -0.1550229 | -0.4353718 | 0.1275793 |
| 4 | 14 | congruent minus blank | -0.0419563 | -0.3971495 | 0.3856576 |
| 5 | 14 | congruent minus blank | -0.1591416 | -0.4995511 | 0.2095518 |
| 6 | 14 | congruent minus blank | -0.1725070 | -0.4610634 | 0.1502083 |
| 1 | 15 | congruent minus blank | -0.1426987 | -0.7984420 | 0.4105439 |
| 2 | 15 | congruent minus blank | -0.3221983 | -0.9105744 | 0.2582093 |
| 3 | 15 | congruent minus blank | 0.0577423 | -0.2971430 | 0.4186432 |
| 4 | 15 | congruent minus blank | 0.0155497 | -0.3833256 | 0.4516939 |
| 5 | 15 | congruent minus blank | 0.0545983 | -0.3267916 | 0.4830101 |
| 6 | 15 | congruent minus blank | 0.0775013 | -0.3923485 | 0.5825811 |
| 1 | 6 | incongruent minus blank | 0.0652574 | 0.0255051 | 0.1160720 |
| 2 | 6 | incongruent minus blank | 0.0333759 | 0.0056439 | 0.0727515 |
| 3 | 6 | incongruent minus blank | 0.0996054 | 0.0518752 | 0.1597846 |
| 4 | 6 | incongruent minus blank | 0.0230215 | 0.0005239 | 0.0571582 |
| 5 | 6 | incongruent minus blank | 0.2522284 | 0.1693830 | 0.3436633 |
| 6 | 6 | incongruent minus blank | 0.0696992 | 0.0300407 | 0.1204411 |
| 1 | 7 | incongruent minus blank | 0.0357948 | -0.0266842 | 0.1038955 |
| 2 | 7 | incongruent minus blank | 0.0000990 | -0.0249165 | 0.0265647 |
| 3 | 7 | incongruent minus blank | 0.1074562 | 0.0415744 | 0.1854842 |
| 4 | 7 | incongruent minus blank | -0.0282598 | -0.0618947 | 0.0066385 |
| 5 | 7 | incongruent minus blank | 0.2406813 | 0.1398935 | 0.3545176 |
| 6 | 7 | incongruent minus blank | 0.1672090 | 0.0938739 | 0.2489809 |
| 1 | 8 | incongruent minus blank | -0.1825103 | -0.2507981 | -0.1167168 |
| 2 | 8 | incongruent minus blank | -0.0710247 | -0.1132772 | -0.0306647 |
| 3 | 8 | incongruent minus blank | -0.1075740 | -0.1626967 | -0.0513064 |
| 4 | 8 | incongruent minus blank | -0.2270295 | -0.2876219 | -0.1701120 |
| 5 | 8 | incongruent minus blank | 0.0055521 | -0.0884583 | 0.1196956 |
| 6 | 8 | incongruent minus blank | 0.0654080 | -0.0049419 | 0.1504645 |
| 1 | 9 | incongruent minus blank | -0.3112378 | -0.4005899 | -0.2185888 |
| 2 | 9 | incongruent minus blank | -0.2736505 | -0.3447802 | -0.1990057 |
| 3 | 9 | incongruent minus blank | -0.2302012 | -0.3088334 | -0.1550389 |
| 4 | 9 | incongruent minus blank | -0.3977470 | -0.4949064 | -0.2949033 |
| 5 | 9 | incongruent minus blank | -0.1493659 | -0.2354242 | -0.0480819 |
| 6 | 9 | incongruent minus blank | -0.0756301 | -0.1386837 | -0.0076847 |
| 1 | 10 | incongruent minus blank | -0.3861038 | -0.5155142 | -0.2513526 |
| 2 | 10 | incongruent minus blank | -0.3649206 | -0.4781384 | -0.2471557 |
| 3 | 10 | incongruent minus blank | -0.3812850 | -0.4854964 | -0.2724611 |
| 4 | 10 | incongruent minus blank | -0.2377836 | -0.3785910 | -0.0896505 |
| 5 | 10 | incongruent minus blank | -0.3155044 | -0.4337350 | -0.1788459 |
| 6 | 10 | incongruent minus blank | -0.2153815 | -0.2883915 | -0.1441249 |
| 1 | 11 | incongruent minus blank | -0.2122111 | -0.4184909 | -0.0097063 |
| 2 | 11 | incongruent minus blank | -0.1983540 | -0.3568964 | -0.0439855 |
| 3 | 11 | incongruent minus blank | -0.2809412 | -0.4151182 | -0.1429701 |
| 4 | 11 | incongruent minus blank | -0.0323623 | -0.2489970 | 0.1882865 |
| 5 | 11 | incongruent minus blank | -0.2772612 | -0.4439921 | -0.0937045 |
| 6 | 11 | incongruent minus blank | -0.2909476 | -0.3900971 | -0.1843613 |
| 1 | 12 | incongruent minus blank | -0.1290275 | -0.3645474 | 0.0962887 |
| 2 | 12 | incongruent minus blank | -0.1370274 | -0.3654997 | 0.0697213 |
| 3 | 12 | incongruent minus blank | -0.1727589 | -0.3331901 | -0.0078045 |
| 4 | 12 | incongruent minus blank | 0.0266732 | -0.2223650 | 0.2964874 |
| 5 | 12 | incongruent minus blank | -0.2490585 | -0.4468342 | -0.0500317 |
| 6 | 12 | incongruent minus blank | -0.2518908 | -0.3809951 | -0.1099974 |
| 1 | 13 | incongruent minus blank | 0.0458617 | -0.2556690 | 0.3635158 |
| 2 | 13 | incongruent minus blank | 0.0085758 | -0.2747837 | 0.3043983 |
| 3 | 13 | incongruent minus blank | -0.0639821 | -0.2784411 | 0.1679976 |
| 4 | 13 | incongruent minus blank | 0.0931941 | -0.2602886 | 0.4403618 |
| 5 | 13 | incongruent minus blank | -0.1461754 | -0.3978298 | 0.1134345 |
| 6 | 13 | incongruent minus blank | -0.4025060 | -0.5734964 | -0.2226895 |
| 1 | 14 | incongruent minus blank | -0.0407115 | -0.3652506 | 0.2661312 |
| 2 | 14 | incongruent minus blank | -0.0049052 | -0.3446538 | 0.3314930 |
| 3 | 14 | incongruent minus blank | -0.1042485 | -0.3485542 | 0.1231077 |
| 4 | 14 | incongruent minus blank | 0.0069196 | -0.3118088 | 0.3377940 |
| 5 | 14 | incongruent minus blank | -0.0969491 | -0.3963804 | 0.1809947 |
| 6 | 14 | incongruent minus blank | -0.1850259 | -0.4240226 | 0.0255308 |
| 1 | 15 | incongruent minus blank | 0.1417184 | -0.5290678 | 0.6361888 |
| 2 | 15 | incongruent minus blank | -0.0637831 | -0.6546504 | 0.4853379 |
| 3 | 15 | incongruent minus blank | 0.0117260 | -0.3003643 | 0.2952493 |
| 4 | 15 | incongruent minus blank | 0.5182057 | -0.0217086 | 0.9386904 |
| 5 | 15 | incongruent minus blank | 0.0774429 | -0.2801920 | 0.4336469 |
| 6 | 15 | incongruent minus blank | -0.2601641 | -0.5911069 | 0.0093791 |
Plot the contrasts in predicted hazard values.
p22 <- ggplot(epreds_pid_diffs, aes(x=timebin, y=diff,
fill=contrast, color=contrast)) +
stat_lineribbon(alpha = 0.4,
point_interval = "mean_qi",
.width=c(.8,.95)) +
#stat_halfeye(point_interval = "mean_qi",.width=c(.8,.95)) +
geom_hline(yintercept = 0, color = "red", lty = 3) +
scale_fill_brewer(palette = "Dark2") +
scale_color_brewer(palette = "Dark2") +
scale_x_discrete(labels = str_c("(",(c(6:15)-1)*40,",",c(6:15)*40,"]"), breaks = 6:15) +
theme(axis.text.x = element_text(angle=90)) +
scale_y_continuous(limits=c(-0.6,0.6)) +
labs(x = "time bin", y = "difference in predicted hazard") +
#ggtitle("Subject-specific AMEs") +
facet_wrap(~pid, ncol=3)
p22## Warning: Removed 1455 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).## Warning: Removed 2308 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).## Warning: Removed 35 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).## Warning: Removed 3578 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).## Warning: Removed 163 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).## Warning: Removed 498 rows containing missing values or values outside the scale range
## (`stat_slabinterval()`).
Combine the plots.
# combine 2 plots
p11_theme <- p11 + theme(legend.position = "none")
(p11_theme / p22) +
plot_annotation(tag_levels = 'A') +
# plot_layout(guides = "collect",
# axes = "collect_x") &
theme(legend.position = "bottom")Save the combined plot.
# save combined plot
ggsave("Tutorial_2_Bayesian/figures/M1i_ame_combined.png", width = 10, height = 12, dpi = 800)Example conclusions for model M1i.
What can we conclude from model M1i about our research question, i.e., the temporal dynamics of the effect of prime-target congruency on RT? In other words, in which of the 40-ms time bins between 200 and 600 ms after target onset does changing the prime from blank to congruent or incongruent affect the hazard of response occurrence (for a prime-target SOA of 187 ms)?
If we want to study the population-level effect of prime type on hazard, uncontaminated by inter-individual differences, we can base our conclusion on Figure 7A. The contrast “congruent minus blank” was estimated to be 0.09 hazard units in bin 6 (95% CrI = [0.02, 0.17]), and 0.14 hazard units in bin 7 (95% CrI = [0.04, 0.25]). For the other bins, the 95% credible interval contained zero. The contrast “incongruent minus blank” was estimated to be 0.09 hazard units in bin 6 (95% CrI = [0.01, 0.21]), -0.19 hazard units in bin 9 (95% CrI = [-0.31, -0.06]), -0.27 hazard units in bin 10 (95% CrI = [-0.45, -0.04]), and -0.23 hazard units in bin 11 (95% CrI = [-0.40, -0.03]). For the other bins, the 95% credible interval contained zero. Note that we could also have calculated hazard ratios instead of hazard differences.
There are thus two phases of performance for the average person between 200 and 600 ms after target onset. In the first phase, the addition of a congruent or incongruent prime stimulus increases the hazard of response occurrence compared to blank prime trials in the time period (200, 240]. In the second phase, only the incongruent prime decreases the hazard of response occurrence compared to blank primes, in the time period (320,440]. The sign of the effect of incongruent primes on the hazard of response occurrrence thus depends on how much waiting time has passed since target onset.
The posterior distribution of each contrast can also be summarized by considering its proportion below or above some value, like zero. For example, here are the proportions that each contrast is larger and smaller than 0:
pabove <- epreds_grand_diffs %>%
group_by(timebin,contrast) %>%
summarize(prop_above = mean(diff > 0)) %>%
arrange(contrast,timebin)## `summarise()` has grouped output by 'timebin'. You can override using the
## `.groups` argument.pbelow <- epreds_grand_diffs %>%
group_by(timebin,contrast) %>%
summarize(prop_below = mean(diff < 0)) %>%
arrange(contrast,timebin)## `summarise()` has grouped output by 'timebin'. You can override using the
## `.groups` argument.table_prop <- pabove %>% inner_join(pbelow, by = c("timebin","contrast"))write_csv(table_prop, file="Tutorial_2_Bayesian/tables/contrasts_props.csv")Show table.
kable(table_prop, caption = "Summarizing the posterior distributions of each contrast by their proportion below and above zero. prop_below = proportion below zero; prop_above = proportion above zero.")## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")## Warning in attr(x, "format"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")| timebin | contrast | prop_above | prop_below |
|---|---|---|---|
| 6 | congruent minus blank | 0.985875 | 0.014125 |
| 7 | congruent minus blank | 0.991000 | 0.009000 |
| 8 | congruent minus blank | 0.951375 | 0.048625 |
| 9 | congruent minus blank | 0.789750 | 0.210250 |
| 10 | congruent minus blank | 0.243000 | 0.757000 |
| 11 | congruent minus blank | 0.154250 | 0.845750 |
| 12 | congruent minus blank | 0.327250 | 0.672750 |
| 13 | congruent minus blank | 0.281125 | 0.718875 |
| 14 | congruent minus blank | 0.193875 | 0.806125 |
| 15 | congruent minus blank | 0.469125 | 0.530875 |
| 6 | incongruent minus blank | 0.983500 | 0.016500 |
| 7 | incongruent minus blank | 0.918250 | 0.081750 |
| 8 | incongruent minus blank | 0.115125 | 0.884875 |
| 9 | incongruent minus blank | 0.003000 | 0.997000 |
| 10 | incongruent minus blank | 0.011125 | 0.988875 |
| 11 | incongruent minus blank | 0.017125 | 0.982875 |
| 12 | incongruent minus blank | 0.040875 | 0.959125 |
| 13 | incongruent minus blank | 0.203375 | 0.796625 |
| 14 | incongruent minus blank | 0.271625 | 0.728375 |
| 15 | incongruent minus blank | 0.578500 | 0.421500 |
Thus, the probability that the contrast “congruent minus blank” is larger than 0, is larger than .9 in bins 6 to 8. And the probability that the contrast “incongruent minus blank” is smaller than 0, is larger than .9 in bins 9 to 12.
If we want to focus more on inter-individual differences, we can study the subject-specific differences in hazard in Figure 7B. Note that three participants (1, 2, and 3) show a negative difference for the contrast “congruent minus incongruent” in bin (360,400] – subject 2 also in bin (320,360].
Future studies could (a) increase the number of participants to estimate the proportion of “dippers” in the subject population, and/or (b) try to explain why this dip occurs. For example, @panisWhatShapingRT2016 concluded that active, top-down, task-guided response inhibition effects emerge around 360 ms after the onset of the stimulus following the prime (here: the target stimulus). Such a top-down inhibitory effect might exist in our priming data set, because after some time participants will learn that the first stimulus is not the one they have to respond to; To prevent a premature overt response to the prime they thus might gradually increase a global response threshold during the remainder of the experiment, which could result in a lower hazard in congruent trials compared to blank trials, for bins after ~360 ms, and towards the end of the experiment. This effect might be masked for incongruent primes by the response competition effect.
Interestingly, all subjects show a tendency in their mean difference (congruent minus blank) to “dip” around that time (Figure 9). Therefore, future modeling efforts could incorporate the trial number into the model formula, in order to also study how the effects of prime type on hazard change on the long experiment-wide time scale, next to the short trial-wide time scale. In Tutorial_2a.Rmd we provide a number of model formula that should get you going.
6. Plot the logit and complementary log-log (cloglog) link functions.
probability <- (1:99999)/100000
logistic <- function(x) { return( 1/(1+exp(-1*x)) )}
logit <- function(x) { return( log(x/(1-x)) )}
inverse_cloglog <- function(x) { return( 1-(exp(-1*exp(x))) )}
cloglog <- function(x) { return( log(-1*log(1-x)) )}
cloglog <- cloglog(probability)
logit <- logit(probability)
dataplot <- cbind(probability,cloglog,logit) %>%
as_tibble() %>%
pivot_longer(cloglog:logit, names_to = "link", values_to = "value")
ggplot() +
geom_hline(yintercept=0, color="grey") +
geom_line(data=dataplot,aes(y=value,x=probability,colour=link),linewidth=1) +
geom_line(data=dataplot,aes(y=value,x=probability,colour=link),linewidth=1) +
scale_color_brewer(palette = "Dark2") +
geom_vline(xintercept = logistic(0), linetype="dotted", linewidth = 0.3) +
geom_vline(xintercept = inverse_cloglog(0), linetype="dotted", linewidth = 0.3) +
annotate("text", x = logistic(0)-.02, y = -6, label = "logistic(0) = 0.5", angle = 90, size=4) +
annotate("text", x = inverse_cloglog(0)+.02, y = -6, label = "inverse_cloglog(0) = 0.6321", angle=90,size=4) +
labs(x = "Probability",
y = "logit or cloglog scale") +
theme(panel.grid = element_blank())
ggsave("Tutorial_2_Bayesian/figures/linkfunctions.png", width = 8, height = 8, dpi = 800)7. Visualize different prior distributions on the logit and cloglog scales.
To gain a sense of what prior logit values would approximate a uniform distribution on the probability (i.e., discrete-time hazard) scale, Solomon Kurz simulated a large number of draws from the Uniform(0,1) distribution, converted those draws to the log-odds metric, and fitted a Student’s t model. Here we do the same for prior cloglog values: simulate a large number of draws from U(0,1), convert them to the cloglog metric, and fit a skew-normal model (due to the asymmetry of the cloglog link function), to gain a sense of what prior cloglog values would approximate a uniform distribution on the probability (i.e., discrete-time hazard) scale.
Simulate, convert, and fit.
set.seed(11)
logit <- function(x) { return( log(x/(1-x)) )}
cloglog <- function(x) { return( log(-1*log(1-x)) )}
# generate draws from U(0,1) and convert
dat <-
tibble(p = runif(1e6, 0, 1)) %>%
mutate(g = logit(p),
c = cloglog(p))
# display
dat %>%
ggplot(aes(x = c)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
# fit model
fit_skewN <-
brm(data = dat,
family = skew_normal(),
c ~ 1,
chains = 4, cores = 4,
file = "Tutorial_2_Bayesian/models/fit_skewN")fit_skewN <- readRDS("Tutorial_2_Bayesian/models/fit_skewN.rds")
summary(fit_skewN) ## Family: skew_normal
## Links: mu = identity; sigma = identity; alpha = identity
## Formula: c ~ 1
## Data: dat (Number of observations: 1000000)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.59 0.00 -0.59 -0.59 1.00 3127 2835
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.26 0.00 1.26 1.26 1.00 2951 2825
## alpha -4.22 0.01 -4.25 -4.19 1.00 1142 1319
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Now we can reverse the process. We simulate from the skew-Normal distribution based on the posterior means for mu, sigma, and alpha, and then convert the results into the probability (i.e., discrete-time hazard) metric.
set.seed(11)
inverse_cloglog <- function(x) { return( 1-(exp(-1*exp(x))) )}
tibble(c = rskew_normal(1e6, mu=-0.59, sigma = 1.26, alpha = -4.22) ) %>%
mutate(p = inverse_cloglog(c)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
Visualize how seven prior distributions on the logit and/or cloglog scales look on the probability scale.
logistic <- function(x) { return( 1/(1+exp(-1*x)) )}
logit <- function(x) { return( log(x/(1-x)) )}
inverse_cloglog <- function(x) { return( 1-(exp(-1*exp(x))) )}
cloglog <- function(x) { return( log(-1*log(1-x)) )}
set.seed(23)
# A N(0,4) prior on the logit and cloglog scales pushes mass to probabilities of 0 and 1
pr1 <- tibble(prior = rnorm(1e6, mean = 0, sd = 4)) %>%
ggplot(aes(x = prior)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
annotate(geom="text", x=-13, y=60000, label="N(0,4)",
color="red", size=4) +
annotate(geom = 'text', label = 'A', x = -Inf, y = Inf, hjust = 0, vjust = 1, size=8)+
theme(panel.grid = element_blank())
l1 <- tibble(log_odds = rnorm(1e6, mean = 0, sd = 4)) %>%
mutate(p = logistic(log_odds)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
ggtitle("logistic(prior)") +
theme(panel.grid = element_blank())
c1 <- tibble(cloglog_prob = rnorm(1e6, mean = 0, sd = 4)) %>%
mutate(p = inverse_cloglog(cloglog_prob)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
ggtitle("inverse_cloglog(prior)")+
theme(panel.grid = element_blank())
# A N(0,2) prior on the logit and cloglog scales pushes mass to probabilities of 0 and/or 1
pr2 <- tibble(prior = rnorm(1e6, mean = 0, sd = 2)) %>%
ggplot(aes(x = prior)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
annotate(geom="text", x=-6, y=60000, label="N(0,2)",
color="red", size=4) +
annotate(geom = 'text', label = 'B', x = -Inf, y = Inf, hjust = 0, vjust = 1, size=8)+
theme(panel.grid = element_blank())
l2 <- tibble(log_odds = rnorm(1e6, mean = 0, sd = 2)) %>%
mutate(p = logistic(log_odds)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
c2 <- tibble(cloglog_prob = rnorm(1e6, mean = 0, sd = 2)) %>%
mutate(p = inverse_cloglog(cloglog_prob)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
# A student-t(df=7.61) prior with scale 1.57 on the logit scale approximates a uniform distribution on the probability scale. This might be a good prior to use for the alpha parameters or Intercept in a logit-hazard model.
pr3 <- tibble(prior = rt(1e6, df = 7.61)* 1.57) %>%
ggplot(aes(x = prior)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
annotate(geom="text", x=-12, y=140000, label="t(7.61,0,1.57)",
color="red", size=4) +
annotate(geom = 'text', label = 'C', x = -Inf, y = Inf, hjust = 0, vjust = 1, size=8)+
theme(panel.grid = element_blank())
l3 <- tibble(log_odds = rt(1e6, df = 7.61)* 1.57) %>%
mutate(p = logistic(log_odds)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
c3 <- tibble(cloglog_prob = rt(1e6, df = 7.61)* 1.57) %>%
mutate(p = inverse_cloglog(cloglog_prob)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
# A Normal(0,1) prior on the logit scale gently regularizes p towards .5.
pr4 <- tibble(prior = rnorm(1e6, mean = 0, sd = 1))%>%
ggplot(aes(x = prior)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
annotate(geom="text", x=-3, y=60000, label="N(0,1)",
color="red", size=4) +
annotate(geom = 'text', label = 'D', x = -Inf, y = Inf, hjust = 0, vjust = 1, size=8)+
theme(panel.grid = element_blank())
l4 <- tibble(log_odds = rnorm(1e6, mean = 0, sd = 1))%>%
mutate(p = logistic(log_odds)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
#geom_vline(xintercept=logistic(0), color="red") +
theme(panel.grid = element_blank())
c4 <- tibble(cloglog_prob = rnorm(1e6, mean = 0, sd = 1))%>%
mutate(p = inverse_cloglog(cloglog_prob)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
# A skew_Normal(-0.59,1.26,-4.22) prior on the cloglog scale approxiates a uniform distr. on the hazard scale. This uninformative prior might be good for the alpha parameters or Intercept in a cloglog-hazard model.
pr5 <- tibble(prior = rskew_normal(1e6, mu=-0.59, sigma = 1.26, alpha = -4.22)) %>%
ggplot(aes(x = prior)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
annotate(geom="text", x=-5.5, y=55000, label="skew_N(-0.59,1.26,-4.22)",
color="red", size=4) +
annotate(geom = 'text', label = 'E', x = -Inf, y = Inf, hjust = 0, vjust = 1, size=8)+
theme(panel.grid = element_blank())
l5 <- ggplot()
c5 <- tibble(cloglog_prob = rskew_normal(1e6, mu=-0.59, sigma = 1.26, alpha = -4.20)) %>%
mutate(p = inverse_cloglog(cloglog_prob)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
# The skew_Normal(-1,1,-2) on the cloglog scale is a weakly informative prior for the alpha parameters or Intercept in a cloglog-hazard model because hazard values below .5 more likely than values above .5 in general.
pr6 <- tibble(prior = rskew_normal(1e6, mu=-1, sigma = 1, alpha = -2)) %>%
ggplot(aes(x = prior)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
annotate(geom="text", x=-5, y=55000, label="skew_N(-1,1,-2)",
color="red", size=4) +
annotate(geom = 'text', label = 'F', x = -Inf, y = Inf, hjust = 0, vjust = 1, size=8)+
theme(panel.grid = element_blank())
l6 <- ggplot()
c6 <- tibble(cloglog_prob = rskew_normal(1e6, mu=-1, sigma = 1, alpha = -2)) %>%
mutate(p = inverse_cloglog(cloglog_prob)) %>%
ggplot(aes(x = p)) +
geom_histogram(bins = 50) +
scale_y_continuous(NULL, breaks = NULL) +
theme(panel.grid = element_blank())
((l1 + pr1 + c1) / (l2 + pr2 + c2) / (l3 + pr3 + c3) / (l4 + pr4 + c4)/ (l5 + pr5 + c5) / (l6 + pr6 + c6) ) &
theme(text = element_text(size = 10, face = "bold"),
title = element_text(size = 10, face = "bold"))
ggsave("Tutorial_2_Bayesian/figures/plot_of_priors.png", width=14, height=13,dpi=800)8. Prior predictive checks
8.1. Model M0i
Following Gelman et al. (2020), we perform a prior predictive check to see if the binary observations generated from the specified prior distribution(s) reflects our prior beliefs. First, sample the prior distributions using sample_prior=“only”.
plan(multicore)
model_M0i_prior <-
brm(data = data_M0i,
family = bernoulli(link="cloglog"),
formula = event ~ 0 + timebin + (0 + timebin | pid),
prior = priors_M0i,
chains = 1, cores = 1,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
sample_prior = "only",
file = "Tutorial_2_Bayesian/models/model_M0i_prior")model_M0i_prior <- readRDS("Tutorial_2_Bayesian/models/model_M0i_prior.rds")Next, use them to predict prior data. But to better understand what the function add_predicted_draws is doing (which we will typically use), we first simulate prior data for one timebin manually.
# prior predictive check for 1 time bin
# see http://bruno.nicenboim.me/bayescogsci/
# cloglog and inverse-cloglog link functions
inverse_cloglog <- function(x) {return(1-(exp(-1*exp(x))))}
cloglog <- function(x) {return(log(-1*log(1-x)))}
# Specify
N_samples <- 2000 # the number of samples from the skew_normal
N_obs <- 100 # the number of simulated observations per sample
set.seed(1)
# function to generate binary observations from samples
cloglog_predictive_distribution <- function(samples, Nobs) {
# empty data frame with headers
df_pred <- tibble(trial = numeric(0),
y_pred = numeric(0),
iter = numeric(0))
# i iterates from 1 to the length of samples
for (i in seq_along(samples)) {
cloglog_haz <- samples[i]
df_pred <- bind_rows(df_pred,
tibble(trial = 1:N_obs,
y_pred = rbinom(Nobs, 1, inverse_cloglog(cloglog_haz)),
iter = i))
}
df_pred
}
# sample prior cloglog-hazard values for 1 bin
cloglog_samples <- rskew_normal(N_samples,-1,1,-2)
# apply the function
prior_pred <- cloglog_predictive_distribution(cloglog_samples,
N_obs)
# calculate hazard as the mean over binary observations
prior_pred_haz <- prior_pred %>%
group_by(iter) %>%
summarise(pred_haz = mean(y_pred))
# plot
ggplot(prior_pred_haz, aes(x=pred_haz)) +
geom_histogram(binwidth=.01) +
scale_x_continuous(limits = c(0,1)) +
labs(x = "simulated hazard",
title = "Prior predictive distribution",
subtitle = str_c(N_samples," samples; ",
N_obs," observations per sample"))ggsave("Tutorial_2_Bayesian/figures/M0i_ppc_1bin.png", width = 10, height = 8, dpi = 800)knitr::include_graphics("Tutorial_2_Bayesian/figures/M0i_ppc_1bin.png")
The shape of the prior predictive distribution for 1 bin looks as expected (compare with Supplementary Figure 2 (row F) in the Supplemental Material): hazard values below .5 are more likely than hazard values above .5.
Now use add_predicted_draws() to simulate prior-based data for our 10 bins and 6 subjects, while taking into account samples from the priors for the standard deviation of the random effects, and correlations between parameters.
# Generate prior predictive hazard functions for 6 participants
newdata_prior = tibble(timebin = 6:15) %>%
expand_grid(pid = 1:6)
# Specify
N_obs = 100
set.seed(2)
# prepare columns of data set
df_pred <- tibble(timebin = integer(0),
pid = integer(0),
.row = integer(0),
.chain = integer(0),
.iteration = integer(0),
.draw = integer(0),
.prediction = integer(0),
obs = integer(0))
# Call add_predicted_draws() for each simulated observation
for(i in 1:N_obs){
prior_pred <- add_predicted_draws(model_M0i_prior,
newdata=newdata_prior,
summary=F,
ndraws=NULL) %>% # all 2000
mutate(obs = i)
df_pred <- bind_rows(df_pred, prior_pred)
}
# calculate hazard per draw (average across observations)
df_pred_haz <- df_pred %>%
group_by(pid,timebin,.draw) %>%
summarise(pred_haz = mean(.prediction))
# plot 200 draws per participant
ggplot(data = df_pred_haz %>% filter(.draw < 200),
aes(x=timebin, y=pred_haz, group=.draw)) +
geom_line(color="black") +
geom_line(data=df_pred_haz %>% filter(.draw == 20),
color="red",
linewidth=2) +
scale_y_continuous(limits = c(0,1)) +
scale_x_continuous(limits = c(6,15), breaks = c(6:15)) +
labs(x = "time bin",
y = "simulated hazard",
title = "Prior predictive distributions",
subtitle = str_c("200 samples; ",N_obs," observations per sample")) +
facet_wrap(~pid)ggsave("Tutorial_2_Bayesian/figures/M0i_ppc.png", width = 10, height = 8, dpi = 800)knitr::include_graphics("Tutorial_2_Bayesian/figures/M0i_ppc.png")
The red hazard function is one of the 200 functions plotted, for each of six participants. They can take on any shape, with values below .5 being more likely than values above .5.
After our prior predictive check, we can fit model M0i to get the posterior distributions (see section 3).
8.2 Model M0c
Perform a prior predictive check to see if the binary observations generated from the specified prior distributions reflect our prior beliefs. To get an idea of how the intercept and both slopes control the shape of the hazard function, the following code allows to change the coefficients and generate (cloglog-)hazard functions. Once you have found the coefficients for the intercept, period_9 and period_9_sq that generate a variety of hazard functions spanning your prior beliefs, you can then use them to determine a suited mean and standard deviation for each coefficient.
dat <- tibble(x = -3:6,
x2 = x^2,
cloglog1 = -1.8 + 0.3*x - 0.1*x2, # black
cloglog2 = -1 + 0.41*x - 0.04*x2, # green
cloglog3 = -3 - 0.13*x + 0.12*x2, # red
cloglog4 = -0.40 + 0.6*x - 0.11*x2, # blue
haz1 = 1 - exp(-1*exp(cloglog1)),
haz2 = 1 - exp(-1*exp(cloglog2)),
haz3 = 1 - exp(-1*exp(cloglog3)),
haz4 = 1 - exp(-1*exp(cloglog4)))
p1<-ggplot(dat, aes(x=x)) +
geom_line(aes(y=cloglog1), color="black") +
geom_line(aes(y=cloglog2), color="green") +
geom_line(aes(y=cloglog3), color="red") +
geom_line(aes(y=cloglog4), color="blue") +
scale_y_continuous(limits=c(-6,2)) +
scale_x_continuous(breaks=c(-3:6)) +
labs(x = "timebin", y = "cloglog-hazard")
p2<-ggplot(dat, aes(x=x))+
geom_line(aes(y=haz1), color="black") +
geom_line(aes(y=haz2), color="green") +
geom_line(aes(y=haz3), color="red") +
geom_line(aes(y=haz4), color="blue") +
scale_y_continuous(limits=c(0,1)) +
scale_x_continuous(breaks=c(-3:6)) +
labs(x = "timebin", y = "hazard")
p1|p2
For example, the coefficients for x vary from -0.13 to 0.6, and we might specify N(.3,.25) as a prior distribution. Those for x-squared vary from -0.11 to 0.12, and we might specify N(0,.06) as a prior distribution.
Suppose we expect (cloglog-)hazard functions that increase almost linearly with time. We also set the standard deviation of the normal prior for class “sd” to a low value, to minimize the effect of the random effects while inspecting the relation between the selected priors for period_9 and period_9_sq and the shape of the prior predictive distributions.
priors_M0c_ppc1 <- c(
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "Intercept"),
set_prior("normal(0.3,.001)", class = "b", coef = "period_9"),
set_prior("normal(0,.001)", class= "b", coef="period_9_sq"),
set_prior("normal(0, .01)", class = "sd"),
set_prior("lkj(2)", class = "cor")
)Sample the prior distributions.
plan(multicore)
model_M0c_prior1 <-
brm(data = data_M0c,
family = bernoulli(link="cloglog"),
formula = event ~ 0 + Intercept + period_9 + period_9_sq +
(0 + Intercept + period_9 + period_9_sq | pid),
prior = priors_M0c_ppc1,
chains = 1, cores = 1,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
sample_prior = "only",
file = "Tutorial_2_Bayesian/models/model_M0c_prior1")model_M0c_prior1 <- readRDS("Tutorial_2_Bayesian/models/model_M0c_prior1.rds")Check the prior distributions.
# extract prior draws
post <- as_draws_df(model_M0c_prior1) %>%
select(-lp__) %>%
as_tibble()
tidy_prior <- post %>%
select(starts_with("b_"), .chain, .iteration, .draw) %>%
rename(chain=.chain, iter=.iteration, draw=.draw) %>%
pivot_longer(-c(chain, draw, iter))
# plot
ggplot(tidy_prior, aes(x = name, y = value, fill = name)) +
stat_halfeye(alpha=0.7,
point_interval = "median_qi",
.width=c(.8,.99),
show.legend = F) +
labs(x = "parameter",
title = "Samples from the prior distributions")
Next, generate prior predictive check for each participant.
Create a function.
plot_ppc <- function(fit_prior, N_samples, N_obs, new_data_prior){
set.seed(1)
# prepare columns of data set with predictions
df_pred <- tibble(period_9 = integer(0),
period_9_sq = integer(0),
pid = integer(0),
.row = integer(0),
.chain = integer(0),
.iteration = integer(0),
.draw = integer(0),
.prediction = integer(0),
obs = integer(0))
# extract predictions
for(i in 1:N_obs){
prior_pred <- add_predicted_draws(fit_prior,
newdata = new_data_prior,
summary=F,
ndraws=N_samples) %>%
mutate(obs = i)
df_pred <- bind_rows(df_pred, prior_pred)
}
# calculate hazard per draw (average across obs)
df_pred_haz <- df_pred %>%
group_by(pid,period_9,.draw) %>%
summarise(pred_haz = mean(.prediction))
# plot N_samples draws per participant, highlighting three
sel = sample(1:N_samples,3,replace=F)
ggplot(df_pred_haz %>%
filter(.draw <N_samples), aes(x=period_9, y=pred_haz, group=.draw)) +
geom_line(color="black") +
# highlight
geom_line(data=df_pred_haz %>% filter(.draw ==sel[1]),
color="red", linewidth=2) +
geom_line(data=df_pred_haz %>% filter(.draw == sel[2]),
color="green", linewidth=2) +
geom_line(data=df_pred_haz %>% filter(.draw == sel[3]),
color="blue", linewidth=2) +
scale_y_continuous(limits = c(0,1)) +
scale_x_continuous(limits = c(-3,6), breaks = c(-3:6)) +
facet_wrap(~pid)
}Apply the function.
# make new data
newdata_prior = tibble(period_9 = c(-3:6)) %>%
mutate(period_9_sq = period_9^2) %>%
expand_grid(pid = 1:6)
# make ppc plot
p1 <- plot_ppc(model_M0c_prior1,2000,100,newdata_prior)
p1For each of 2000 draws from the prior distributions, a black hazard function is plotted. Three randomly selected hazard functions are highlighted in red, green, and blue.
ggsave("Tutorial_2_Bayesian/figures/M0c_ppc1.png", width = 10, height = 8, dpi = 800)knitr::include_graphics("Tutorial_2_Bayesian/figures/M0c_ppc1.png")
Or suppose we believe apriori that the hazard functions are flat. We specify our new priors:
priors_M0c_ppc2 <- c(
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "Intercept"),
set_prior("normal(0,.001)", class = "b", coef = "period_9"),
set_prior("normal(0,.001)", class= "b", coef="period_9_sq"),
set_prior("normal(0, .001)", class = "sd"),
set_prior("lkj(2)", class = "cor")
)Sample the prior distributions 2000 times:
plan(multicore)
model_M0c_prior2 <-
brm(data = data_M0c,
family = bernoulli(link="cloglog"),
formula = event ~ 0 + Intercept + period_9 + period_9_sq +
(0 + Intercept + period_9 + period_9_sq | pid),
prior = priors_M0c_ppc2,
chains = 1, cores = 1,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
sample_prior = "only",
file = "Tutorial_2_Bayesian/models/model_M0c_prior2")model_M0c_prior2 <- readRDS("Tutorial_2_Bayesian/models/model_M0c_prior2.rds")Apply the function plot_ppc().
# make new data
newdata_prior = tibble(period_9 = c(-3:6)) %>%
mutate(period_9_sq = period_9^2) %>%
expand_grid(pid = 1:6)
# make ppc plot
p2 <- plot_ppc(model_M0c_prior2,2000,100,newdata_prior)
p2Now the prior predictive hazard functions are all rather flat, with an intercept controlled by the prior skew-normal distribution, i.e., values below .5 are more likely than values above .5.
ggsave("Tutorial_2_Bayesian/figures/M0c_ppc2.png", width = 10, height = 8, dpi = 800)knitr::include_graphics("Tutorial_2_Bayesian/figures/M0c_ppc2.png")
Finally, perform a ppc for the selected prior distributions for M0c (see lines 220-224), and increase the standard deviation for the normal prior for class “sd” somewhat.
priors_M0c_ppc3 <- c(
set_prior("skew_normal(-1,1,-2)", class = "b", coef = "Intercept"),
set_prior("normal(0.3,.25)", class = "b", coef = "period_9"),
set_prior("normal(0,.06)", class= "b", coef="period_9_sq"),
set_prior("normal(0, .1)", class = "sd"),
set_prior("lkj(2)", class = "cor")
)plan(multicore)
model_M0c_prior3 <-
brm(data = data_M0c,
family = bernoulli(link="cloglog"),
formula = event ~ 0 + Intercept + period_9 + period_9_sq +
(0 + Intercept + period_9 + period_9_sq | pid),
prior = priors_M0c_ppc3,
chains = 1, cores = 1,
iter = 3000, warmup = 1000,
control = list(adapt_delta = 0.999,
step_size = 0.04,
max_treedepth = 12),
seed = 12, init = "0",
sample_prior = "only",
file = "Tutorial_2_Bayesian/models/model_M0c_prior3")model_M0c_prior3 <- readRDS("Tutorial_2_Bayesian/models/model_M0c_prior3.rds")Apply the function plot_ppc().
# make new data
newdata_prior = tibble(period_9 = c(-3:6)) %>%
mutate(period_9_sq = period_9^2) %>%
expand_grid(pid = 1:6)
# make ppc plot
p3 <- plot_ppc(model_M0c_prior3,2000,100,newdata_prior)
p3A wide range of differently shaped hazard functions can be observed, and early bins are unlikely to have high hazard values, consistent with our prior beliefs. Note that some simulated hazard functions start off with a high hazard in bin -3. These reflect a combination of negative slopes and a positive coefficient for period_9_sq. These possibilities are required to encompass hazard functions that initially decrease and then increase, as observed in the empirical hazard functions for incongruent primes.
ggsave("Tutorial_2_Bayesian/figures/M0c_ppc3.png", width = 10, height = 8, dpi = 800)knitr::include_graphics("Tutorial_2_Bayesian/figures/M0c_ppc3.png")